libname in "INSERT PATH TO SAS LIBRARY";
options linesize = 72;

*
   This run first carries out the MDS of the 13
   political figures from the 2004 election, using
   the LOS dissimilarities. Then, the SAS macro,
   "%boot" creates 50 bootstrap replications of
   the LOS matrix, carries out the MDS on each one,
   and then uses the Schonemann and Carroll
   procedure to rotate each bootstrap replication
   of the MDS configuration to optimal fit 
   with the "real" MDS configuration.
;

data los;
input name $8. gbush kerry nader 
cheney edwards lbush hclint bclint 
powell ashcroft mccain dempty reppty;
cards;
gwbush         0.0 73.0 62  8.0 68.0 20.0 51.5 41.0 24.0  7 25.5 50  5.0
kerry         73.0  0.0 56 78.0  1.0 54.0 15.0 17.0 47.0 77 37.0  2 74.5
nader         62.0 56.0  0 72.0 59.0 53.0 60.0 49.0 58.0 70 39.0 57 71.0
cheney         8.0 78.0 72  0.0 74.5 25.5 65.0 51.5 29.0 12 30.0 66  4.0
edwards       68.0  1.0 59 74.5  0.0 44.0 14.0 16.0 46.0 76 38.0  3 69.0
lbush         20.0 54.0 53 25.5 44.0  0.0 42.0 34.0  9.5 23 22.0 45 18.0
hclinton      51.5 15.0 60 65.0 14.0 42.0  0.0 19.0 32.0 67 40.0 13 55.0
bclinton      41.0 17.0 49 51.5 16.0 34.0 19.0  0.0 31.0 61 36.0 11 48.0
powell        24.0 47.0 58 29.0 46.0  9.5 32.0 31.0  0.0 28  9.5 35 21.0
ashcroft       7.0 77.0 70 12.0 76.0 23.0 67.0 61.0 28.0  0 33.0 63  6.0
mccain        25.5 37.0 39 30.0 38.0 22.0 40.0 36.0  9.5 33  0.0 43 27.0
dempty        50.0  2.0 57 66.0  3.0 45.0 13.0 11.0 35.0 63 43.0  0 64.0
reppty         5.0 74.5 71  4.0 69.0 18.0 55.0 48.0 21.0  6 27.0 64  0.0
;
run;

proc mds 
   coef = identity   condition = matrix
   data = los       dimension = 2
   fit = distance    formula = 1
   out = mdsc
   level = ordinal   untie
   pfinal;
id name;
var gbush kerry nader 
   cheney edwards lbush hclint bclint 
   powell ashcroft mccain dempty reppty;
title "MDS, creating target configuration";
run;

data mdsconf;
set mdsc;
if _type_ = "CONFIG";
keep _name_ dim1 dim2;
run;

proc print data = mdsconf;
run;

data therms;
set in.therms;
array flip gwbush--reppty;
   do over flip;
   flip = 100 - flip;
   end;
run;

proc means data = therms;
title "  ";
run;

%macro boot (nbs);

%do booti = 1 %to &nbs;

proc iml worksize = 300;

start readdata;
use therms;
read all var _num_ into therms04 [colname = names];

value = (10 + &booti) # (&booti # 5);
index = ceil(ranuni(j(nrow(therms04),1,value)) # nrow(therms04) );
bdata = therms04[index,];

create bdata from bdata [colname = names];
append from bdata;
finish;

run readdata;
quit;

OPTIONS PAGESIZE = 58;
PROC IML WORKSIZE=300;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          This PROC IML routine calculates the line-of-            */
/*          sight measure of interobject similarity.  The            */
/*          reference for this procedure is:                         */
/*                                                                   */
/*                Rabinowitz, George (1976) A Procedure for          */
/*                Ordering Object Pairs Consistent with the          */
/*                Multidimensional Unfolding Model.                  */
/*                PSYCHOMETRIKA 41: 349-373.                         */
/*                                                                   */
/*          The proper citation for using this SAS procedure is:     */
/*                                                                   */
/*                Jacoby, William G. (1993) A SAS/IML Macro          */
/*                for Calculating the Line-of-Sight Measure          */
/*                of Interobject Dissimilarity.                      */
/*                PSYCHOMETRIKA 58: 511-512.                         */
/*                                                                   */
/*          This procedure assumes an input SAS data set that        */
/*          consists of N subjects' ratings of K objects.  Hence     */
/*          The input data are an N by K matrix.  There should       */
/*          be no missing values, and the scores should be           */
/*          set so that SMALLER values are more positive ratings,    */
/*          and HIGHER values are less positive ratings.             */
/*                                                                   */
/*-------------------------------------------------------------------*/

START NAME;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Module NAME obtains the names of the variables           */
/*          from the input data set, by reading in a single          */
/*          observation, using the COLNAME specification, and        */
/*          then freeing the data.  The names of the rating          */
/*          variables are then stored in the row vector, NAMES.      */
/*          By default, the routine uses the most recently           */
/*          created SAS dataset for input data.  A different         */
/*          SAS dataset can be specified in the USE statement.       */
/*                                                                   */
/*-------------------------------------------------------------------*/

USE bdata;
READ POINT 1 VAR _NUM_ INTO DUMMY [COLNAME=NAMES];
FREE DUMMY;
FINISH;

START READIN;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Module READIN reads the data from the input SAS          */
/*          dataset, and calculates the pairwise sums and            */
/*          differences from the rating scores.  Begin by            */
/*          looping through the columns of the input dataset.        */
/*                                                                   */
/*                I1:  Index for first var in current pair           */
/*             NAME1:  Name of first var in current pair             */
/*               ONE:  Vector of ratings of first object             */
/*                     in the current pair                           */       
/*                                                                   */
/*-------------------------------------------------------------------*/

   DO I1 = 1 TO NCOL(NAMES) - 1;
   NAME1 = NAMES[,I1];
   READ ALL VAR NAME1 INTO ONE;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          I2:  Index for second var in current pair                */
/*       NAME2:  Name of second var in current pair                  */
/*                                                                   */
/*-------------------------------------------------------------------*/

     DO I2 = I1 + 1 TO NCOL(NAMES);
     NAME2 = NAMES[,I2];
     READ ALL VAR NAME2 INTO TWO;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Create vector of current sums (SUMVEC) and               */
/*          vector of current differences (DIFVEC).  Then            */
/*          sort SUMVEC from smallest to largest, and sort           */
/*          DIFVEC from largest to smallest.                         */
/*                                                                   */
/*-------------------------------------------------------------------*/

     SUMVEC = ONE + TWO;
     DIFVEC = ABS(ONE - TWO);
     FREE TWO;
     SORTVEC = SUMVEC;
     SUMVEC[RANK(SUMVEC),] = SORTVEC;
     SORTVEC = DIFVEC;
     DIFVEC[(NROW(DIFVEC) + 1) - RANK(DIFVEC),] = SORTVEC;
     FREE SORTVEC;
       IF I1 = 1 & I2 = 2 THEN DO;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          At first pair, create matrix of cumulative, sorted       */
/*          pairwise sums and differences.  NOBS is the total        */
/*          number of observations in the input data set, and        */
/*          NROWS is the maximum number of iterations that will      */
/*          be used in this run.  After the current sums and         */
/*          differences are concatenated to the SUM and DIF          */
/*          matrices, SUMVEC and DIFVEC can be freed.  At the        */
/*          end of the READIN module, both SUM and DIF will have     */
/*          NROWS rows and K(K-1)/2 columns.                         */
/*                                                                   */
/*-------------------------------------------------------------------*/

       NOBS = NROW(SUMVEC);
       NROWS = ROUND (NOBS / 25);
       SUM = CUSUM(SUMVEC[1:NROWS,]);
       DIF = CUSUM(DIFVEC[1:NROWS,]);
       FREE SUMVEC DIFVEC;
       END;
       IF I1 > 1 | I2 > 2 THEN DO;
       SUM = SUM || CUSUM(SUMVEC[1:NROWS,]);
       DIF = DIF || CUSUM(DIFVEC[1:NROWS,]);
       FREE SUMVEC DIFVEC;
       END;
     END;
   FREE ONE;
   END;
FINISH;

START CALC;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Module CALC uses the cumulative, sorted SUM and          */
/*          DIF matrices to calculate the LOS matrix.                */
/*          Matrices are defined as follows:                         */
/*                                                                   */
/*          CURSTATS:  Vector of current LOS statistics              */
/*          STATNAME:  Vector of labels for printed output           */
/*            NOINCR:  Number of iterations needed to check          */
/*                     for no increase in ADEQU value                */
/*                                                                   */
/*-------------------------------------------------------------------*/       


   CURSTATS = J(1,3,0);
   STATNAME = {DISCRIM DENSE ADEQU};
   NOINCR = ROUND(NOBS / 150);
   IF NOINCR < 4 THEN NOINCR = 4;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Start looping through rows of pairwise sums and          */
/*          differences.  Sum corresponding entries together,        */
/*          obtain rank order of resultant sums.  Then use           */
/*          DISNUM and NDIST to obtain number of distinct            */
/*          ranks on this iteration.                                 */
/*                                                                   */
/*-------------------------------------------------------------------*/

   DO I = 1 TO NROWS;
   DISNUM = RANKTIE(SUM[I,] + DIF[I,]);
   NDIST = 0;
     DO J = 1 TO NCOL(DISNUM);
     IF DISNUM[1,J] > 0 THEN NDIST = NDIST + 1;
       IF J < NCOL(DISNUM) THEN DO;
         DO J2 = J+1 TO NCOL(DISNUM);
         IF DISNUM[1,J2] = DISNUM[1,J] THEN
            DISNUM[1,J2] = -1;
         END;
       END;
     END;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          CURSTATS contains the LOS statistics for the             */
/*          current row of the SUM and DIF matrices.  These          */
/*          values are concatenated onto the bottom of the           */
/*          STATS matrix, for later printing the iteration           */
/*          history.  CURMAX is the current maximum value of         */
/*          the ADEQU statistic.  The row in which this max          */
/*          ADEQU value occurs is indexed by the var CURLOC.         */
/*          The LOS statistics for the current max ADEQU value       */
/*          are contained in the vector MAXSTATS.  On each           */
/*          iteration, check the ADEQU for that row, and replace     */
/*          max ADEQU information, if necessary.                     */ 
/*                                                                   */
/*-------------------------------------------------------------------*/   

   CURSTATS[,1] = NDIST / (NCOL(DISNUM) - 1);
   CURSTATS[,2] = (NOBS - I) / (NOBS - 1);
   CURSTATS[,3] = CURSTATS[,1] # (CURSTATS[,2] ## 3);
     IF I = 1 THEN DO;
     STATS = CURSTATS;
     CURMAX = CURSTATS[,3];
     CURLOC = I;
     MAXSTATS = CURSTATS;
     END;
     IF I > 1 THEN STATS = STATS // CURSTATS;
     IF CURSTATS[,3] > CURMAX THEN DO;
     CURMAX = CURSTATS[,3];
     CURLOC = I;
     MAXSTATS=CURSTATS;
     END;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Check to see whether ADEQU has not increased over        */
/*          the last NOINCR rows.  If it has not increased, then     */
/*          terminate the procedure.                                 */
/*                                                                   */
/*-------------------------------------------------------------------*/

     IF I >= NOINCR THEN DO;
     C1 = STATS[(I-NOINCR+1):(I-1),3];
     C2 = STATS[(I-NOINCR+2):(I),3];
     IF SUM((C1 - C2) < 0) = 0 THEN GO TO TERM2;
     END;
   END;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          TERM1 is used if the procedure goes through all          */
/*          NROWS rows of the matrices.  The I value is              */
/*          decreased by 1, so it prints the correct number          */
/*          of iterations.                                           */
/*                                                                   */
/*-------------------------------------------------------------------*/

TERM1:
   I = I - 1;
   PRINT  'FINAL ITERATION NUMBER:' I,
     'ROUTINE TERMINATED BECAUSE IT REACHED MAX N OF ITERATIONS,,';
   GO TO DONE;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          TERM2 is used if the ADEQU value converges.              */
/*                                                                   */
/*-------------------------------------------------------------------*/

TERM2:
   PRINT  'FINAL ITERATION NUMBER:' I,
     'ROUTINE TERMINATED BECAUSE ADEQUACY VALUE CONVERGED,,';
   GO TO DONE;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*        DONE obtains the final ranking, based on combined sums     */
/*        and differences, and places it into the row vector,        */
/*        OBJRANK.  It also calculates the Spearman correlation      */
/*        between the ranking based on sums only, and that based     */
/*        on differences only (RHO).  It also prints out the         */
/*        iteration history and the LOS statistics for the           */
/*        final solution.                                            */
/*                                                                   */
/*-------------------------------------------------------------------*/

DONE:
   OBJRANK = RANKTIE(SUM[CURLOC,] + DIF[CURLOC,]);
   RHO = 1 - ((6 # SSQ(RANKTIE(SUM[CURLOC,]) - RANKTIE(DIF[CURLOC,])))
         /(NCOL(SUM) # ((NCOL(SUM)##2) - 1)));
   ITER = (1:NROW(STATS))`;   
   ITNAME = {'ITER #'};
   PRINT 'ITERATION HISTORY', ITER [COLNAME=ITNAME FORMAT=3.] 
     STATS [COLNAME = STATNAME FORMAT=7.3],,
     'FINAL LOS STATISTICS:',,
     'SELECTED OBJECT RANKING FROM ITERATION #:' CURLOC ,,
     MAXSTATS [COLNAME=STATNAME FORMAT=7.3] ,,
     'SUM-DIFF RHO:' RHO [FORMAT=7.3];
   FINISH;

START SHAPE;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Module SHAPE places the final LOS values, obtained       */
/*          from the OBJRANK vector, into a square                   */
/*          K by K matrix called DISSIM.                             */
/*                                                                   */
/*-------------------------------------------------------------------*/
                
DISSIM = J(NCOL(NAMES),NCOL(NAMES),0);
I = 0;
   DO I1 = 1 TO NCOL(NAMES) - 1;
     DO I2 = I1 + 1 TO NCOL(NAMES);
     I = I + 1;
     DISSIM[I1,I2] = OBJRANK[,I];
     DISSIM[I2,I1] = OBJRANK[,I];
     END;
   END;
PRINT 'FINAL MATRIX OF DISSIMILARITIES',,
      DISSIM [ROWNAME=NAMES COLNAME=NAMES];
FINISH;

START WRITE;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Module WRITE creates a new SAS dataset, called           */
/*          DISSIM, from the previously-obtained DISSIM matrix.      */
/*          This dataset (with K observations and K variables)       */
/*          contains the pairwise LOS dissimilarity values.  The     */
/*          variable names are identical to those in the input       */
/*          SAS dataset.  Also, remember that the LOS values are     */
/*          dissimilarities-- lower values are more similar pairs    */
/*          and higher values are more dissimilar pairs              */
/*                                                                   */
/*-------------------------------------------------------------------*/

CREATE DISSIM FROM DISSIM [COLNAME = NAMES];
APPEND FROM DISSIM;
FINISH;

RUN NAME;
RUN READIN;
RUN CALC;
RUN SHAPE;
RUN WRITE;

proc mds 
   coef = identity   condition = matrix
   data = dissim     dimension = 2
   fit = distance    formula = 1
   initial = mdsc
   noprint
   out = ebs
   level = ordinal   untie
   pfinal;
var gwbush kerry nader 
   cheney edwards lbush hclinton bclinton 
   powell ashcroft mccain dempty reppty;
title "MDS of bootstrap sample";
run;

data estims;
set ebs;
if _type_ = "CONFIG";
keep _name_ dim1 dim2;
run;

proc iml worksize = 300;
*
   Module readin reads the 
   target configuration and
   the bootstrap replication
   configuration, which will
   be rotated to optimal fit with
   the target.
;
start readin;

use mdsconf;
read all var _num_ into target 
   [colname = dims];

use estims;
read all var _num_ into bconf
   [colname = dims];

read all var _char_ into name;

*
   corr1 is the correlation between
   the target matrix and the bootstrap
   replication matrix, before 
   rotation
;
corr1 = (sum(bconf # target)) /
   sqrt(ssq(bconf) # ssq(target));

finish;

*
   Module rotate will rotate the
   bootstrap replication matrix to
   optimal fit with the target
;
start rotate;

cmatrix = target` * bconf;

call svd(pmat, phimat, qmat, cmatrix);

tmatrix = qmat * pmat`;

bconf2 = bconf * tmatrix;

finish;

*
   Module dilate will dilate
   the bootstrap replication
   matrix to optimal conformity
   with the target
;
start dilate;
numer = trace(target` * bconf2);
denom = trace(bconf` * bconf);

sval = numer / denom;

bconf3 = bconf2 # sval;

finish;
*
   Module translate shifts
   the origin in the rotated,
   dilated space to come into
   maximum conformity with the
   target
;
start translat;
tvec = (nrow(target)##-1) #
   ((target - bconf3)` * j(nrow(target), 1, 1));

bconf4 = bconf3 + (j(nrow(target), 1, 1) * tvec`);

corr4 = (sum(bconf4 # target)) /
   sqrt(ssq(bconf4) # ssq(target));

finish;

start writerep;
correl = corr1 || corr4;
corname = {"corr1" "corr4"};

create correl from correl [colname = corname];
append from correl;

create estim from bconf4 [rowname = name colname = dims];
append from bconf4 [rowname = name];
finish;

run readin;
run rotate;
run dilate;
run translat;
run writerep;

data bsrep;
   do i = 1 to 13;
   replic = &booti;
   output;
   end;
keep replic;
run;

data current;
merge estim bsrep;
run;

   %if &booti = 1 %then %do;
   data resample;
   set current;
   run;

   data correls;
   set correl;
   run;
   %end;

   %else %do;
   data resample;
   set resample current;
   run;

   data correls;
   set correls correl;
   run;
   %end;

%end;
%mend boot;

%boot(50)

proc print data = resample;
title "print of final coordinate replic dataset";
run;

proc print data = correls;
title "Print of final congruence coeff dataset";
run;

